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Monte Carlo simulations of the field substellar mass function (MF) are presented, based 
on the latest brown dwarf evolutionary models from Burrows et al. (1997) and Baraffe et al. 
(2003). Starting from various representations of the MF below 0.1 M Q and the stellar birth rate, 
luminosity functions (LFs) and T e ff distributions are produced for comparison with observed 
samples. These distributions exhibit distinct minima in the mid-type L dwarf regime followed 
by a rise in number density for fainter/cooler brown dwarfs, predicting many more T-type and 
cooler brown dwarfs in the field even for relatively shallow mass functions. Deuterium-burning 
brown dwarfs (0.012 M Q < M < 0.075 M ) dominate field objects with 400 < T e// < 2000 K, 
while non- fusing brown dwarfs make up a substantial proportion of field dwarfs with T e ff < 
500 K. The shape of the substellar LF is fairly consistent for various assumptions of the Galactic 
birth rate, choice of evolutionary model, and adopted age and mass ranges, particularly for field T 
dwarfs, which as a population provide the best constraints for the field substellar MF. Exceptions 
include a depletion of objects with 1200 < T e ff < 2000 K in "halo" systems (ages > 9 Gyr), 
and a substantial increase in the number of very cool brown dwarfs for lower minimum formation 
masses. Unresolved multiple systems tend to enhance features in the observed LF and may 
contribute significantly to the space density of very cool brown dwarfs. However, these effects are 
small (< 10% for T e // < 300 K) for binary fractions typical for brown dwarf systems (10-20%). 
An analytic approximation to correct the observed space density for unresolved multiple systems 
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in a magnitude-limited survey is derived. As an exercise, surface densities as a function of T e f f 
are computed for shallow near- infrared (e.g., 2MASS) and deep red-optical (e.g., UDF) surveys 
based on the simulated LFs and empirical absolute magnitude/T e // relations. These calculations 
indicate that a handful of L and T dwarfs, as well as late-type M and L halo subdwarfs, should 
be present in the UDF field depending on the underlying MF and disk scale height. These 
simulations and their dependencies on various factors provide a means for extracting the field 
substellar MF from observed samples, an issue pursued using 2MASS T dwarf discoveries in 
Paper II. 

Subject headings: Galaxy: stellar content — methods: numerical — stars: low mass, brown 
dwarfs — stars: luminosity function, mass function 
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1. Introduction 

The stellar Initial Mass Function (IMF) is a fundamental quantity in astrophysics. Defined as the total 
number density of stars ever created in a particular environment per unit mass (Miller & Scalo 1979), the 
IMF is a sensitive probe of the star formation process, accounts for the mass budget and evolution of galaxies, 
and determines the evolution of chemical abundances over time. Pioneering work by Salpeter (1955) showed 
that the IMF for field stars in the Solar Neighborhood (masses 0.4 < M < 10 M Q ) could be adequately 
reproduced by a power law, \I/(M) = oc M -2 - 35 , a result that generally persists to this day (Scalo 1998; 
Kroupa 2001; Reid, Gizis, & Hawley 2002). Since that time, many studies of the IMF have been undertaken 
for low and high mass stars of differing populations, and in various regions of the Galaxy and external star 
clusters. Excellent reviews can be found in Miller & Scalo (1979); Scalo (1986); Kroupa (1998); Scalo (1998); 
Reid & Hawley (2000); and Chabricr (2003). 

The IMF is a particularly key measurement in the study of brown dwarfs. These objects comprise the 
low-mass tail of the stellar population, but differ in that they lack sufficient mass to sustain core Hydrogen 
fusion (Hayashi & Nakano 1963; Kumar 1963). Because of this, brown dwarfs never reach the Hydrogen 
main sequence, but instead continually evolve to cooler temperatures and fainter magnitudes. The intrinsic 
faintness of brown dwarfs made them an early candidate for dark matter (Tarter 1975; Bahcall 1984); indeed, 
an extrapolation of the Salpeter IMF yields nearly twice as much mass in brown dwarfs (0.005 < M < 0.075 
M Q ) as in stars (0.075 < M < 40 M Q ). However, number counts of field M dwarfs show a flattening in 
the IMF around 0.3-0.5 M Q (Sandage 1957; Schmidt 1959; Miller & Scalo 1979), and it is now quite clear 
that brown dwarfs are not prolific enough to be the constituents of dark matter. Nevertheless, the number 
density of brown dwarfs may still be a significant fraction or multiple of the stellar density (Reid et al. 1999; 
Kroupa 2001; Chabrier 2002), and the nearest systems to the Sun may in fact be unidentified substellar 
ones. Furthermore, quantifying the IMF in the substellar regime enables a unique exploration of the star 
formation process; in particular, its efficiency at small masses and the lower limit at which self-gravitating 
"stars" can form. 

The IMF is not an observable quantity, and is generally derived from the Luminosity Function (LF), 
$(Mb ;), the number density of stars observed in a defined region per unit luminosity. The LF is converted 
into the Present Day Mass Function (PDMF, the number density of stars currently present in a defined 
region per unit mass), using empirical (e.g., Henry & McCarthy 1993) or theoretical (e.g., Baraffe et al. 
1998) mass-luminosity (M-L) relations. For the lowest-mass stars and brown dwarfs (M < 0.1 M©) in well- 
defined regions of space, and assuming no evolution of the star-forming process over time (although its rate 
can change), the IMF is identical to the PDMF and can be referred to simply as the Mass Function (MF). 

While this technique is suitable for low-mass stellar populations, substellar MF determinations are hin- 
dered by their thermal evolution. A brown dwarf with an observed luminosity and/or effective temperature 
(T e ff) has a wide range of possible masses depending on its age. This mass-age degeneracy is not critical 
for young cluster brown dwarf populations, where members are assumed to be approximately coeval (e.g., 
White & Ghez 2001). In the Galactic disk, however, stars and brown dwarfs can span a fairly broad range of 
ages, from a few tens of Myr to <~ 10 Gyr. In other words, there is no single M-L relation that can be used to 
convert the LF into the MF for brown dwarfs in the field. Field brown dwarfs are also generally older than 
their young cluster counterparts, so that the lowest mass field objects can be exceedingly faint, requiring 
deep and/or wide area surveys to detect sufficient numbers. Nevertheless, the physical properties of evolved 
brown dwarfs are better understood than their younger counterparts, without the complications of youth- 
ful accretion or rapid evolution. Furthermore, the nearby population of stars is not affected by reddening; 



- 3 - 



can be more easily followed-up with spectroscopic, parallactic, and high-resolution imaging observations (to 
derive physical characteristics and multiplicity); and, assuming that it is well-mixed, is generally devoid of 
foreground or background contamination. 

This article is the first of a two-part series investigating the substellar MF in the Solar Neighborhood, 
by comparing simulated LFs to a magnitude-limited sample of T dwarfs (Burgasser et al. 2003a) identified 
in the Two Micron All Sky Survey (Cutri ct al. 2003, hereafter 2MASS). T dwarfs are a spectroscopic 
class of brown dwarfs that exhibit CH4 absorption (Burgasser et al. 2002b; Geballe et al. 2002), implying 
T e ff < 1300 K (Kirkpatrick et al. 2000; Golimowski et al. 2004). In this article, Monte Carlo simulations 
of the field substellar MF are examined, and dependencies on various input parameters are investigated. 
These simulations are comparable to those of Allen et al. (2004), who constrain the substellar MF through 
Bayesian techniques. The implementation of the simulations presented here is described in § 2, which 
includes discussion of the various input distributions and evolutionary models used. An in-depth analysis 
of the derived LF and T e // distributions and their features is given in § 3. In § 4, the sensitivity of these 
distributions between the evolutionary models employed, different birth rates, different age and mass limits, 
and the influence of unresolved multiple systems is explored. Surface density predictions based on the 
simulations are derived for both shallow and deep magnitude-limited surveys in § 5. Results are summarized 
in § 6. 



2. The Simulations 

2.1. General Description of the Problem 

The purpose of these simulations is to create a statistical link between the MF and LF, or more gen- 
erally a link between the fundamental properties of brown dwarfs - mass, age, and metallicity - and their 
observables - T e // and luminosity. This link is made through evolutionary models coupled to non-grey 
model atmospheres. In this study, we assume that all brown dwarfs are described by a single distribution 
for each of their fundamental parameters, denoted P (x), where x is the fundamental property in question 1 . 
The fundamental distributions examined are summarized in Table 1 and described in detail below. 

Following traditional practice, it is assumed that the MF does not evolve with time, so that the number 
density of stars ever created per unit mass and per unit time, C(M,t) (termed the Creation Function by 
Miller & Scalo 1979) can be separated into mass- and time-dependent functions: 

C(M,t) = tf(M)x6(t)/T (1) 

(Miller & Scalo 1979), where ^(M) is the MF, b(t) is the birth rate (number density of stars born per unit 
time), and T a is the age of the Galaxy, assumed here to be 10 Gyr. The separation of C(M,t) enables the 
examination of the MF and birth rates separately. The metallicity (Z) and mass ratio (q) distributions are 
also assumed to be independent so that they may be treated separately as well. Clearly, these assumptions 
may not accurately reflect the detailed stellar formation history of the Galaxy (e.g., chemical evolution is 
ignored; see Edvardsson ct al. 1993), but are adequate for current substellar MF determinations. 



x Note that the mass distribution is the MF; i.e., P(M) = *(M). 
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2.2. Fundamental Distributions 

2.2.1. The Mass Distribution 
Six MFs were examined, including five power-law distributions: 

*(M) oc M~ a (2) 
with a = 0.0, 0.5, 1.0, 1.5, and 2.0; and the lognormal MF from Chabricr (2001): 

(log M-log M c ) 2 

*(logM) oc e ^ (3) 

(see also Miller & Scalo 1979), where M c = 0.1 M Q and a = 0.627. Note that *(logM) = In 10 M f(M). 
The baseline simulations incorporate a mass range 0.01 < M < 0.1 M Q , with the upper limit set by the 
evolutionary models and the lower limit set to provide enough objects in the higher mass bins, particularly 
for the steeper power laws. Lower mass limits ranging from 0.001 to 0.015 M Q were also examined in order 
to measure the influence of a minimum "cut-off" mass (M mi „) in brown dwarf formation (§ 4.2.4). 



2.2.2. The Age Distribution 



The age distribution is related to the birth rate 2 by P(t) = b(T a — t). The birth rate does not influence 
the shape of the MF if the latter is assumed not to evolve. However, as brown dwarfs themselves evolve 
thermally, the age distribution can influence the LF. Five birth rates were examined, as illustrated in Figure 
1: 
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The first ("constant") birth rate is the simplest and most frequently employed for MF simulations, and a 
number of authors have asserted its legitimacy based on studies of the Galactic star formation history (SFH) . 
Miller & Scalo (1979) argue that the SFH must be roughly flat over the age of the Galaxy to explain the 



2 Herc, t denotes the age of an object, counting backwards from the current epoch. The birth rate 6(r) is generally defined 
in terms of increasing time r = T — t. 
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continuity of the MF between low and high mass stars; formation rates, kinematics, and spatial distribution 
of planetary nebulae and white dwarfs; nucleosynthesis yields; distribution of H II regions; and theoretical 
predictions at that time. Soderblom, Duncan, & Johnson (1991) claim no evidence of variation of the star 
formation rate over the past 10 9 years based on the activity distribution of G and K stars, (although a 
reanalysis by Rocha-Pinto & Maciel (1998) argues otherwise; see below). Boissier & Prantzos (1999) also 
find little evidence for variation between recent and early SFHs based on the metallicity distribution of 
G-dwarfs. 

The second ("exponential") birth rate has been used to model Galactic star formation (Tinsley 1974; 
Miller & Scalo 1979) because of its simple form. This birthrate is consistent with a star formation rate that 
scales with the average gas density (Miller & Scalo 1979), with an e-folding time r g — 5 Myr; i.e., half of the 
age of the Galaxy (as adopted here). Miller & Scalo (1979) find this rate to be marginally consistent with 
continuity arguments, although in general there are no empirical data that strongly support this function. 

More recent studies have suggested that the SFH is not strictly monotonic, but can be characterized 
by a series of burst events. Barry (1988) point out an apparent increase in star formation 400 Myr ago, a 
result supported by an examination of the white dwarf luminosity function by Noh & Scalo (1990, however, 
see Soderblom, Duncan, & Johnson 1991). The presence of perhaps three burst episodes in the SFH of the 
Galaxy is detailed in Majewski (1993). To model such a non-monotonic birth rate, a smoothed version of the 
empirical results of Rocha-Pinto et al. (2000) was used, based on the chromospheric ages of late- type dwarfs. 
This "empirical" birth rate exhibits peaks 0—1, 2-5, and 7-9 Gyr ago, with somewhat lower formation rates 
in between these bursts. As discussed in Rocha-Pinto et al. (2000), a smoothed distribution may hide more 
dramatic swings in the Galactic SFH, but the detection of such events are hindered by uncertainties in stellar 
age measurements. Note that this non-monotonic birth rate does not violate the continuity arguments of 
Miller & Scalo (1979). 

The fourth birth rate examined is a novel one assuming star formation has occurred entirely in young 
clusters, in a series of short-lived formation bursts evenly and randomly distributed over the age of the 
Galaxy. The formation period is short in young clusters, t c ; w 10 — 20 Myr (White & Ghez 2001), so that 
this birth rate approximates a "stochastic" formation process. A total of N c i = 50 clusters was assumed, 
randomly distributed over the age of the Galaxy and each producing an equal number of brown dwarfs over 
the same formation timescale. The birth rate distribution of each cluster was assumed to be Gaussian with 
a characteristic time scale r c ; = 10 Myr. These assumptions arc not necessarily representative of the true 
yields and lifetimes of young clusters in the Galaxy, but are suitable for this study. 

Finally, the fifth ("halo") birth rate considers only brown dwarfs born within a 1 Gyr burst 9 Gyr in 
the past, and is meant to represent the conditions of the Galactic halo or old globular cluster substcllar 
populations (Reid & Hawley 2000). 

For each of these birth rates, an age range of 0.01 < t < 10 Gyr is nominally adopted, although minimum 
ages of 1 to 100 Myr were also examined to investigate the contribution of young populations in the simulated 
LFs (see § 4.2.3). 

2.2.3. The Metallicity Distribution 

The choice of a metallicity distribution is primarily constrained by the evolutionary models used (§ 
2.3), both of which assume solar abundances. Therefore, a constant distribution P{Z) = 1 is adopted with 
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Z = Zq. This choice is supported by the fact that 70% of disk stars have abundances —0.3 < [m/H] < 0.15 
(Reid & Hawley 2000), but requires that there is no significant contamination by other Galactic populations 
(e.g., thick disk and halo brown dwarfs) in the observed sample. 



To convert our fundamental properties to observables, we used the most recent evolutionary calculations 
from the Tucson (Burrows et al. 1997) and Lyon (Baraffe et al. 2003) groups. Both of these models employ 
non-grey atmospheres in which condensate opacity is ignored (so-called "COND" models; Allard et al. 2001), 
largely consistent with the observed spectra of mid-type M and mid- and late-type T dwarfs (Tsuji, Ohnaka, 
& Aoki 1996). Chabrier et al. (2000a) have also derived evolutionary tracks for "DUSTY" atmosphere 
models, which retain condensate material in their atmosphere, more appropriate for warmer late-type M 
and L dwarfs (Tsuji, Ohnaka, & Aoki 1996). However, these authors find <10% difference in the evolution 
of luminosity and T e // between the COND and DUSTY models. This is a relatively small deviation given 
the potentially larger systematic uncertainties arising from the complex evolution of condensates in cool M 
and L dwarf atmospheres (Ackerman & Marley 2001; Burgasser et al. 2002a; Tsuji 2002; Cooper et al. 2003) 
and current observational uncertainties (Burgasser 2001; Cruz et al. 2003). DUSTY evolutionary tracks are 
therefore ignored in this investigation. 

In Figure 1, the evolution of T e // with time for the two sets of models employed are compared for 
masses 0.001 < M < 0.1 M Q and ages 1 Myr to 10 Gyr. Over much of this parameter space evolutionary 
tracks are consistent to within 10%, with the Baraffe models predicting slightly higher temperatures at a 
particular mass and age for M < 0.06 M© and lower temperatures for M > 0.08 M©. At early ages (t < 5 
Myr) the Burrows models are significantly hotter (20-25%) for M > 0.06 M Q . At later ages (t > 5 Gyr), the 
two models again deviate significantly (20-35%) for 0.06 < M < 0.08 M©, with the Baraffe models being 
both hotter and more luminous. This is due to the higher Hydrogen Burning Minimum Mass (HBMM) for 
the Burrows models, 0.075 versus 0.072 M©. Finally, the Burrows tracks diverge more substantially around 
the HBMM, with a difference of 1500 K between 0.075 and 0.09 M at 10 Gyr, as compared to 600 K for 
identical masses and ages in the Baraffe models. The effects of these differences on the simulated LFs are 
described in § 4.2.1. 



For each simulation, a total of 3xl0 6 objects were assigned a set of fundamental properties x% (x = M, 
t) by selecting the random parameter & from a uniform distribution over the range:, 



2.3. Evolutionary Models 



2.4. Implementation of the Simulations 



& e [min(P{x}),max(P{a;})], 
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Each object was then assigned a luminosity (and hence Mj, ; = —2.5 log (L/Lq) +4.74; Livingston 2000) and 
T e ff by linearly interpolating the grid points of the appropriate evolutionary model to that object's mass 
and logarithmic age. Minimum values of Mboi = 50 mag and T e ff = 100 K were assigned to the derived 
observables if the model values fell below these limits. The observable distributions $(Mb ;) (in units of 
pc~ 3 mag -1 ) and <&(T e //) (in units of pc~ 3 [100 K] -1 ) were then determined by binning the observable 
parameters every 0.5 mag and 100 K, respectively. 

In order to extract meaningful comparisons between the various distributions and empirical data, sim- 
ulated MF number densities were normalized to the mean of the field low-mass star (0.1-1.0 M ) MFs 
of Reid et al. (1999) and Chabricr (2001), *(M) = 0.35( 57 ^)" 113 and *(M) = 0.67( 57 ^)- 1 - 55 pc -3 
Mq 1 , respectively. Over the range 0.09-0.1 M Q , these mass functions yield an average number density of 
0.0055±0.0018 pc~ 3 . The number of objects over the same mass range in each simulation sample was nor- 
malized to this value, and that normalization applied to each output distribution. The 30% discrepancy 
between the two stellar mass functions at 0.1 M is significant, but as all of the distributions are scaled 
by this factor, adjustment to refined estimates of the low-mass stellar space density can be readily made. 
Values for W(M) for each of the MFs employed are given in Table 2. 

Figure 3 shows the resulting MF distributions for simulations with baseline parameters: Baraffe evo- 
lutionary models, P(t) — constant, 0.01 < M < 0.1 M Q , and 0.01 < t < 10 Gyr. These distributions are 
consistent with their analytic forms to better than 3% over most of the mass range examined, with somewhat 
larger scatter (not exceeding 10%) in the lowest mass bins for the steepest power law distributions. These 
accuracies are identical for simulations using lower cutoff masses. Therefore, numerical uncertainties are neg- 
ligible in comparison to, e.g., observational uncertainties (Paper II) and differences between the evolutionary 
models (§4.2.1). 

3. Results 

A total of 32 Monte Carlo simulations were run to examine the various input parameters described above. 
Resulting observable distributions for the baseline simulations are given in Tables 3 and 4 and diagrammed 
in Figures 4 and 5. 

3.1. The Luminosity Function 

Figure 4 diagrams the derived LFs, Q(Mboi), for baseline parameters and for each of the MFs examined. 
Also labelled in this plot (and subsequent figures) are the approximate Mtois for spectral types (SpT) L0, 
L5, T5, and T8, based on empirical measurements by Golimowski et al. (2004). At bright magnitudes, there 
is a peak at M\, i <~ 13 (SpT < L0) which is less pronounced for the steeper MFs but yields the same density 
of objects (0.01 pc -3 mag -1 ) for all MFs. This peak is almost entirely comprised of low mass stars (0.08 
< M < 0.1 M Q ), and the fixed density reflects the adopted normalization. The drop off in $(M;, ;) toward 
brighter luminosities is an artifact of the upper mass limit (0.1 M ) of the simulations. At M^oi ~ 15 (SpT 
~ L5) there is a local minimum in 5>(Mf, ;), a feature that has also been seen in the simulations of Chabricr 
(2003) and Allen et al. (2004, their "Trough B"). The origin of this trough may be seen in the divergence of 
the evolutionary tracks in Figure 2 around T e ff ~ 1800 K (corresponding to M^oi ~ 15). At late ages, this 
temperature straddles the HBMM, and hence most brown dwarfs have cooled to lower temperatures and 
fainter luminosities. Sources older than 1 Gyr tend to dominate the overall population for a flat birthrate 
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(see § 4.1); hence, the narrow range of masses sampling these luminosities at late ages implies fewer sources 
overall. Note that shallower power laws produce a more pronounced trough. Toward fainter magnitudes, 
&(Mb i) rises, more significantly for steeper power laws due to the greater proportion of low mass (and hence 
intrinsically fainter for a given age) brown dwarfs. Each distribution exhibits a broad peak at these faint 
magnitudes, with the location of the maximum depending on the steepness of the MF: Mboi <~ 18 for a 
= and Mboi ~ 22 for a = 2. Indeed, beyond Mboi ~ 18 (SpT <~ T7), there is a substantial increase in 
the contrast between the various MFs, with up to 25 times more brown dwarfs between a = 2 and at 
Mboi ~ 21. The lognormal $(Mboi) lies between those of the a — 0.5 and 1.0 MFs, and is generally flat 
between 18 < Mboi < 21. Below Mboi ~ 22, there is a steep drop off in all of the $(Mf, j) distributions due to 
both the adopted lower mass limit (0.01 M Q for the simulations diagrammed in Figure 4) and the adopted 
maximum age (10 Gyr). A lower minimum mass and/or and an older population would result in a turnover 
in the LF at fainter magnitudes. 



3.2. The T eff Distribution 

Figure 5 compares the <&(T e jj) distributions for the same simulations. The bright magnitude peak seen 
in the $(M& j) distribution is evident at T e ff <~ 2500— 2700 K, although it likely underestimates the actual 
number of stars/brown dwarfs at these higher temperatures because of the 0.1 M Q upper mass cutoff. The 
trough in $(Mf, ;) is also seen here, albeit somewhat less pronounced, around 1800-2000 K (SpT <~ L3-L5), 
again due to the rapid cooling of brown dwarfs at these temperatures. At lower T e ffS, all of the distributions 
rise, with the steeper power laws yielding at least an order of magnitude more cold brown dwarfs (T e / / ~ 500 
K) than warm ones (T e// ~ 2000 K). At 1000 K (SpT ~ T6), there is a factor of 8 difference between a = 
and 2, and a factor of 30 difference at 500 K. The resulting densities of cold brown dwarfs are fairly high, 
predicting roughly 25 brown dwarfs with 400 < T^eff ^$ 800 K within 5 pc of the Sun for a = 1. This is 
somewhat less than half the number of main sequence stars in an equivalent volume (Reid, Gizis, & Hawley 
2002). Below 300 K, there is a sharp turnover in $(T e //) similar to that seen in <J>(M ho /) for M bol < 22. 

4. Analysis 

4.1. Composition of $(M bo i) and *(T e// ) 

It is instructive to break down the luminosity and T e ff distributions by mass and age in order to 
examine in detail the origins of the various features seen. Figure 6 shows <J>(T e //) for the a = 0.5 simulation 
for which a low mass cutoff of 0.001 M Q was used (see § 4.2.3). This distribution is broken down into 
groupings of low mass stars (0.075 < M < 0.1 M Q ), Deuterium-burning brown dwarfs (0.012 < M < 0.075 
M ), and non-fusing brown dwarfs (0.001 < M < 0.012 M©). It is clear that the high temperature peak 
in the LF is indeed dominated by main sequence low-mass stars down to 1900-2000 K (SpT ~ L3), with 
a smaller contribution of predominantly young Deuterium-burning brown dwarfs. At cooler temperatures, 
Deuterium-burning brown dwarfs are the dominant population down to T e ff ~ 500 K, encompassing all 
of the currently known field brown dwarfs. Non-fusing brown dwarfs only make a significant contribution 
below this temperature. This segregation of masses in the <E>(T e //) distribution is seen for all of the MFs 
examined. 

An alternate way to examine the mass composition of $(Mb i) and $(T e //) is by computing the median 
mass per luminosity or T e // bin, as diagrammed in Figure 7 for simulations with M min = 0.001 M© and a = 
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0.5 and 1.5. The most likely range of masses in each bin was chosen to comprise 63% of all objects about the 
median value, equivalent to ±ler in a Gaussian distribution. Three trends are immediately discernable; first, 
the median mass decreases toward lower luminosities and cooler temperatures, consistent with the fact that 
lower-mass brown dwarfs start off cooler, and therefore remain cooler, than their higher-mass counterparts 
at any given age. Second, as the median mass relations cross the HBMM, they diverge for different MFs, 
with the steeper distributions exhibiting lower median masses at a given luminosity or temperature. This 
is simply due to the larger number of lower-mass brown dwarfs in the steeper MFs contributing to each of 
the luminosity and temperature bins. Finally, there is a wide range of masses that comprise each luminosity 
and temperature bin, a range that increases for steeper MFs with the inclusion of more low-mass sources. 
In one sense, these substantial mass "uncertainties" highlights the motivation for the simulations — the 
non-unique nature of the field substcllar M-L relation — and demonstrates the substantial uncertainty in 
assigning masses to field objects without age information. 

Figure 8 plots the median age as a function of luminosity and T e ff for the same MF simulations; the 
indicated typical range of ages was computed as above. In this case, the spread in ages in each bin is 
substantial; it is not possible to assign a statistical age with uncertainty better than a few Gyr based on 
luminosity and T e // alone. However, there are some subtle trends in these relations that may have statistical 
merit. There is a notable drop in the median age at the same locations as the troughs in the Q(M}, i) and 
§(T e ff) distributions, around 1500 < T e ff < 2000 K. These features are related, as the higher luminosities 
and hence more rapid evolution of brown dwarfs at these temperatures implies both fewer objects present at 
any given time and very few brown dwarfs remaining or reaching these temperatures at later ages. At earlier 
times, this temperature region encompasses a much broader range of masses and hence a larger percentage 
of the young population. Allen et al. (2004) note a similar age bias amongst L dwarfs in their simulations. 
One consequence of this feature is that L dwarfs in the field should be younger on average than T dwarfs. 
There is some empirical evidence of this form tangential velocity measurements (Vrba et al. 2004) and the 
mass-age-activity trends of late- type M and L dwarfs (Gizis et al. 2000). However, it is important to stress 
that the range of ages sampled at these temperatures is still very large, and individual age determinations 
cannot be precisely determined. The apparent decrease in median age for steeper power-law MFs is again 
due to the greater contribution of lower-mass brown dwarfs, which appear in the higher temperature and 
luminosity bins when they are younger and less evolved. 



4.2. Variations in $ Distributions from Various Factors 

The observable distributions presented above are based primarily on the baseline parameters of a flat 
birth rate, 0.01 < M < 0.1 M Q , 0.01 < t < 10 Gyr, and the Baraffe et al. (2003) evolutionary models. The 
trends identified in these distributions indicate methods of constraining the substellar MF by comparison to 
empirical data; however, they may be confused by other details such as the choice of evolutionary model, the 
form of the birth rate, the age range of field brown dwarfs, and the minimum formation mass. Quantifying 
the influence of these parameters on the shape and scale of the observable distributions provides a measure 
of the systematic uncertainty in the derived MF when comparing to empirical data. 
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4-2.1. Variations due to Choice of Evolutionary Model 

Figure 9 compares <J>(M hoi ) and $(T e //) between the Burrows et al. (1997) and Baraffe et al. (2003) 
models for the a = 0.5 and 1.5 power-law MFs. Some of the variations in the tracks seen in Figure 2, 
particularly at the stellar/substcllar boundary, are reflected in the resulting observable distributions. Most 
notably, the bright peak at Mf> ; ~ 13 {T e ff ~ 2500 K) is far less pronounced in the Burrows et al. (1997) 
model simulations. Low mass stars are instead piled up at slightly brighter luminosities (Mf, ; ~ 12.5) and 
hotter temperatures (T e ff ~ 2800 K). Furthermore, the Burrows et al. (1997) models predict fewer objects 
overall at brighter magnitudes (Mboi < 15) and hotter temperatures (T e ff ^ 1500 K), and more objects 
at fainter magnitudes (M^oi > 22) and colder temperatures (T e ff 400 K) than the Baraffe et al. (2003) 
models. In the T dwarf regime, however, the two sets of models are in fairly good agreement for all of the 
MFs examined. Therefore, the choice of evolutionary model does not appear to affect the interpretation of 
the T dwarf field LF, but can be an important source of systematic uncertainty when examining the LF of 
hotter (L- and M-type) brown dwarfs. 

4-2.2. Variations due to Differing Birth Rates 

Figure 10 compares $(Mf, ;) and (fr(T e ff) for the five birthrates for an a — 0.5 baseline simulation. For 
both distributions, there is effectively no difference between the constant, empirical, and stochastic birth 
rates, the more realistic realizations for the Galactic field population. This result is consistent with the 
findings of Allen et al. (2004), who discern minimal variations in derived LFs between birth rates that are 
constant, based on field star ages (Soderblom, Duncan, & Johnson 1991), and based on star formation rates 
as a function of redshift (Pascual et al. 2001). Therefore, the underlying birth rate generally has a negligible 
effect on the determination of the field MF. 

The more extreme exponential and halo birth rates, however, do modulate the observable distributions, 
with a far more pronounced dip at M boi - 15 - 17 (T e// - 1200 - 2000 K; SpT - T5-L3) and more 
fainter/cooler brown dwarfs. Both of these effects are due to the larger proportion of older, and therefore 
more evolved and fainter, brown dwarfs produced by these birth rates. The differences are most pronounced 
for the halo age distribution, which predicts a substantial deficiency of T e ff ~ 1200 — 2000 K brown dwarfs, 
comprised primarily of mid- and late-type L and early T dwarfs. Note that this deficiency is likely to be 
more pronounced in a real halo population, as the reduced mctallicity typical for halo dwarfs (Gizis 1997) 
imply more transparent atmospheres, enhanced luminosities, and hence more rapid cooling (Burrows et al. 
2001). 

It is interesting to note that all of the distributions are generally consistent between 18 < M^oi < 21 
(500 < T e ff < 1000 K), which encompasses mid- and late-type T dwarfs. These consistencies suggest that 
while the field T dwarf population may be highly sensitive to the underlying MF (§ 3.1), it is generally 
insensitive to the birth rate. In contrast, the L dwarf field population is somewhat less sensitive to the MF 
but may be an excellent probe of extreme Galactic birth rates. These trends are also seen in the steeper 
MFs. 
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4-2.3. Variations due to Differing Age Limits 

Figure 11 compares $(M;, ;) and <&(T e //) for minimum ages of 1 to 100 Myr for the a = 0.5 and 1.5 
MFs. No significant differences are seen between these LFs, primarily because very young (t < 100 Myr) 
brown dwarfs contribute minimally (~1% for a flat birthrate) to the 10 Gyr field LF over the mass range 
examined. Hence, the minimum age of the substellar field population only influences the observed LF if it 
is of order 1 Gyr or later, as seen with the halo birthrate discussed above. 

4-2.4- Variations due to Differing Minimum Mass Cutoffs 

One of the key parameters for low-mass star formation theories is the minimum formation mass, which 
depends not only on the thermodynamical conditions of the initial gas reservoir (through the Jean's mass, 
Jeans 1902), but also on the efficiency and history of accretion early in the formation process. Figure 12 
compares $>(Mi, i) and <&(T e //) for minimum formation masses M m i„ = 0.001, 0.010, and 0.015 M© and the 
a = 0.5 and 1.5 MFs. As expected, reducing M m i„ results in many more intrinsically faint objects, and the 
low-temperature turnover in <E>(T e //) (Figure 5) is essentially absent for M min = 0.001 M Q . However, the 
differences between these distributions are negligible for M bo [ < 20 and Teff > 500 K for both power-law 
MFs. This is consistent with the mass breakdown of <J>(T e ^y) in Figure 6, which shows that the lowest mass 
brown dwarfs contribute significantly only to the lowest temperature/luminosity bins. Thus, the signature 
of a minimum brown dwarf formation mass, unless it is larger than 0.015 M , cannot be detected in the 
currently known sample of field brown dwarfs, which extend only to T e // <~ 700 K (Golimowski et al. 2004; 
Vrba et al. 2004). Determining M mi „ from field measurements will require the discovery of substantially 
cooler brown dwarfs. 



4.3. The Influence of Multiplicity 

Any observed sample of stars or brown dwarfs may include some percentage of unresolved multiple 
systems. Indeed, stellar multiples are more frequent amongst Solar-mass stars than single systems (Abt & 
Levy 1976; Duquennoy & Mayor 1991, ~ 60%), and this frequency may be even higher during the early 
T-Tauri phase (Ghez, Neugebauer, & Mattews 1993). Recent high-resolution imaging studies of low-mass 
stars and brown dwarfs have shown that a small fraction of these systems (^10-20%) are closely-separated 
(a < 15 AU) binaries (Kocrncr et al. 1999; Rcid et al. 2001; Close et al. 2002; Bouy et al. 2003; Burgasser et 
al. 2003b; Gizis et al. 2003), and at least one substellar spectroscopic binary is also known (Basri & Martin 
1999a). All of these systems are unresolved in wide-field imaging surveys such as 2MASS and SDSS; hence, 
brown dwarf samples drawn from those surveys tend to measure the systemic LF, & sys , rather than the 
distribution of individuals, Qi n d- The latter is a more appropriate constraint for star formation theory or 
the Galactic mass budget. 

Unresolved multiple systems induce two effects on the LF: (1) an increase in the number of individ- 
ual sources in the sample; and (2) an increase in the effective volume sampled for each of the individual 
components of a multiple system, as their unresolved combined light allows them to be detected to greater 
distances 3 . To examine these effects, a second series of Monte Carlo simulations were performed, building 



3 An alternate interpretation of this second effect, discussed with the author by L N. Reid (2004, priv. comm.), is that the 
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from the <j>(T e //) distributions from the MF simulations. It was assumed that the sample under investigation 
is a magnitude-limited one for which T e //s can be measured (e.g., through spectral typing or color), the 
typical scenario for samples large enough to measure the mass function. Furthermore, only coeval binary 
systems were considered, and it was assumed that the binary fraction higher order systems are ignored) 

and mass ratio (q = M2/M1) distribution (P(q)) are fixed and independent of mass, and can therefore be 
treated separately from \P(M) and P(t). Finally, it was assumed that the primary of each system has the 
same temperature (T^) as the observed systemic T e ff, as would be the case if the unresolved spectrum 
(and hence spectral type) is dominated by the brighter component. 

An analytic approximation to this problem, appropriate for the overall space density of a population, is 
given in the Appendix. For the simulations, the effect of multiplicity on the observed $(T e //) = <& S ys(T e //) 
was considered by determining the correction factor $j„d(T e //)/$ sys (T e //). First, $ ai/s (T e //) was scaled 
by (1 — fbin) to give the number density of single objects. The remaining fraction of binary systems were 
modelled using N = 10 6 test sources per T e ff bin. Mass ratios for the binaries were assigned from three 
choices of P{q), listed in Table 1 and diagrammed in Figure 13, where q was allowed to vary from 0.001 to 
1. The first ("flat") distribution is generally consistent with results from closely-separated (spectroscopic) 
stellar binary studies (e.g., Mazeh et al. 1992, 2003). The second ("exponential") distribution assumes a 
greater percentage of equal-mass systems, a form consistent with recent studies of low-mass star and brown 
dwarf binaries (Reid et al. 2001; Gizis et al. 2003; Goldberg, Mazeh, & Latham 2003). A value of q c = 0.26 
was derived from a fit to the apparent q distribution of known L and T dwarf binaries (Reid et al. 2001; 
Bouy et al. 2003; Burgasser et al. 2003b; Gizis et al. 2003, Figure 13). Note that this empirical distribution 
has not been corrected for selection effects (e.g., incompleteness for low-g systems), and is therefore purely 
an exploratory one. The third distribution assumes both primaries and secondaries are drawn from the 
same underlying MF, an interpretation put forth by Duquennoy & Mayor (1991) to explain the mass ratio 
distribution of G and K stars (see also Kroupa & Burkcrt 2001). The distribution shown in Figure 13, 
which peaks at lower mass ratios, was generated from Monte Carlo simulations of 10 6 primaries and 10 6 
secondaries, both drawn from an a = 0.5 power-law mass function with M m i n = 0.001 M Q , and random 
pairing (q is fixed to be no greater than unity). 

For each binary simulation, effective temperatures for the secondary components of the binaries, T);^, 
were determined as: 

T ( 2 ) /t(2)\ 1 / 4 

^f-Q < 13 > 

L eff V 7 

where it is assumed that the primary and secondary radii are equal (true to within 10-15% for ages greater 
than 1 Gyr) and L oc M 2 64 (Burrows et al. 2001). Equation 13 matches theoretical evolutionary models 
fairly well but can overestimate T^ by 10-20% for systems that staddle the H- or D-burning limits. It is, 
however, a useful analytical approximation. The N sec = fu n N secondaries were then binned by their T e //, 
and the space density of both primaries and secondaries added to the single star $(T e //) distribution after 
scaling by the factor: 

n(T i= T e// ) 

e(T e// ) = -J- £ (l + qf- 6 T 3/2 , (14) 

-''sec 

z 

where the sum is over the n simulated primaries or secondaries for which Tj = T e ff (after binning). The 



spectrophotomctric distance for an unresolved binary is underestimated due to that system's brighter combined light, resulting 
in an overestimated space density. If the sample is constrained to be volume-limited, the correction for this effect is identical 
to that derived here. 
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factor (1 + qf - 64 )- 3 / 2 = €i compensates for the increased volume sampled by the combined light of the binary 
system (see Appendix, Eqn. A3). Thus, for each T e ff bin k, the temperature distribution of individual 
sources is 

*ind(Tfc) = (l-/wn)*,»a(Tfc) 

+ /wn*«»«(Tfc)e(T fc ) (15) 
+ fbinEj *sy S (T< 1J |Tf = T fc )e(Tf ) 

(2) 

where the final term incorporates secondaries for which —Tfc, but normalizes to the systemic space 
density corresponding to the temperature of the primary. 

Figure 14 plots ^md/^sys as a function of T e ff for various assumptions of fb in , P(q), and a. All of 
distributions show significant structure, with correction values less than unity for Teff iS 2300 K (SpT <~ 
L0), large correction values for T e// ~ 1300-1700 K (SpT ~ T0-L5), a dip around T eff ~ 900 K (SpT ~ 
T7), and a rapid rise toward the coolest temperatures. Most of this structure can be attributed to a trickle- 
down effect amongst the binary secondaries. The low values of <&i n d/$sys for hotter dwarfs is due to the 
lack of secondaries from systems with primaries having M > 0.1 M0 in these simulations; hence, only equal- 
mass/equal-magnitudc binaries contribute. As discussed in the Appendix, a mass-ratio distribution skewed 
toward q = 1 causes the increased effective volumes of these systems to be a larger effect than the addition of 
secondaries into the sample. This feature, then, is an artifact of our simulation upper mass limit. The peak 
at ^eff ~ 1300-1700 K is enhanced by both the decline of the underlying LF at these temperatures and the 
addition of low- and moderate-g secondaries associated with primaries from the 2500-2700 K $(T e ff) peak. 
The valley in $>i n d/$sys a t cooler temperatures is the result of the opposite trend: a rise in the underlying 
LF at these temperatures coincident with a paucity of hotter host primaries hosting low- and moderate-g 
secondaries. At the coolest temperatures, very low-mass secondaries associated with primaries across the 
LF contribute to &i n d 7 and a large correction factor is needed to account for these systems. Note that 
mass ratio distributions skewed toward higher mass ratios (e.g., the exponential P(q)) result in a smaller 
correction factor beyond T e ff ~ 1800 K, due to the paucity of low-mass secondaries contributing to the 
lower temperature bins. A similar trend is seen for steeper MFs, but in this case is due to the steeper rise 
of the underlying LF at lower temperatures. 

The morphology of &ind/®sys implies that unresolved multiplicity tends to enhance key features in the 
LF, resulting in a larger contrast between L and T dwarf numbers in the observed & sys . On the other hand, 
the very coolest and faintest brown dwarfs in an observed sample are largely hidden as low-g secondaries, 
resulting in an artificial flattening of the LF. It is important to note, however, that these variations are 
generally small. For binary fractions typical for brown dwarf systems, |1 — <&i n d/^sys\ < 10% for T e ff > 300 
K, similar to the systematic uncertainties from the evolutionary models and far more accurate than current 
LF measurements of late- type field dwarfs (e.g., Cruz et al. 2003). Even binary fractions as high as 50% 
cause only a 20% shift in the LF in the late L dwarf regime. Hence, the influence of unresolved multiplicity 
will be difficult to discern in the field LF given the current precision of observations; however, as larger field 
samples are generated, the features described above could provide an independent means of probing the mass 
ratio distribution of low mass stars and brown dwarfs. 



5. Surface Density Predictions 

The purpose of the simulations presented here is to place constraints on the substellar MF using LF 
measurements of field brown dwarfs, an issue that will be pursued in detail in Paper II. The simulations can 
alternately be used as a predictive tool; specifically, as a means of estimating the number of brown dwarfs 
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detectable in a particular field imaging survey. To illustrate this, surface densities (E) as a function of T e ff 
for two types of imaging surveys were examined: a shallow near-infrared survey similar to 2MASS, and a 
deep red-optical survey at high Galactic latitude, similar to the Great Observatories Origins Deep Survey 
(Giavalisco 2004) or the Hubble Ultra Deep Field (UDF, Beckwith et al. 2003). 

Starting from the 3>(T e //) distributions with M min = 0.001 M Q , surface densities for a shallow survey 
were computed by assuming that the space density is constant throughout the volume observed, so that 

£afc(T e //) - l(^) 2 ^eff)d 3 max (T eff ). (16) 

Here, d max (T e f j) = iQ- - 2 ( M ( T ^/f)~ m nm)+ 1 p C j s th e limiting detection distance for brown dwarfs in the 
survey to an apparent magnitude limit mu m , and M(T e //) is the absolute magnitude/T e // relation for the 
imaging filter used. The latter can be determined using either theoretical models or empirical data. For a 
deep, high Galactic latitude survey, d max can be comparable to the scale height of the Galactic disk (H z ), so 
that the vertical distribution of sources must be considered. At a height z above/below the Galactic plane, 
the space density of stars scales as 

$(z) = <i> sech 2 (i^-) (17) 

(Rcid & Hawlcy 2000), where <I> is the local space density. The resulting surface density can therefore be 
written as S<jeep = l^sh, where 

MT ,„ } . i^-'^^m^, (18) 

amax(Leff) Jo lti z 

is the scale height correction factor for a survey field at Galactic latitude j3. Note that corrections to the 
radial distribution of sources, which is important for deep surveys extending to several kpc scales close in 
the Galactic plane, are not considered here 4 . 

For the shallow imaging case, ii-band 5 imaging down to mu m = 16 is considered, similar to the sen- 
sitivity limits of 2MASS; while a z' field to mu m — 28 (AB magnitudes) at (3 = 54?5 is considered for the 
deep imaging case, appropriate for the Hubble UDF. My and M z > versus T e // relations for single M6-T8 
dwarfs were determined empirically using photometry compiled from Dahn et al. (2002) and Knapp et al. 
(2004); parallax measurements from Dahn et al. (2002), Tinncy, Burgasser, & Kirkpatrick (2003), and Vrba 
ct al. (2004); and T e // determinations from Golimowski et al. (2004). Linear fits of absolute photometry 
versus log 10 T eff yield: 

M H (T eff ) = 50.27- 11.65 log 10 T e// (a = 0.13 mags) (19) 

for 42 sources with 700 < T eff < 2900 K, and 

M z ,(T eff ) = 45.60 -8.96 log 10 T e// (a = 0.20 mags) (20) 

for 17 sources with 900 < T e ff < 1750 K. These linear relations are extrapolated over the full T e // range 
of our sample. 



4 Intcrstcllar dust absorption would also have a profound effect on Galactic plane surveys, perhaps more so than the radial 
limits of the disk (I. N. Rcid 2004, priv. comm.). The effect of interstellar absorption perpendicular to the plane is ignored here. 

5 The J-band/T e ^ relation for L and T dwarfs exhibits non-monotonic behavior at the L/T transition (Dahn et al. 2002), 
possibly due to the evolution of dust clouds at these temperatures (Burgasser et al. 2002a). This makes the correction from 
T e f f to Mj degenerate. The somewhat less sensitive //-band is therefore used for this exercise. 
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Figure 15 plots S(T e //) for both of the surveys considered. For the shallow case, three populations were 
examined: "disk" dwarfs with mass functions a = 0.5 and 1.5; and a halo birthrate population with a = 1.5 
scaled by a factor of 0.3%, consistent with the relative number of halo to disk stars in the Solar Neighborhood 
(Reid & Hawlcy 2000). In all cases, £ decreases rapidly with T e ff, and the relative densities for the two 
disk MFs scale with the underlying LF. Halo stars arc greatly outnumbered by disk stars, consistent with 
the adopted normalization, but this contrast increases to ~1% in the T dwarf regime. The low densities in 
the L and T dwarf regime (<~ 10 -5 — 10~ 3 deg -2 [100 K] _1 ) imply that substantial areas must be imaged 
to identify a statistically significant number of sources. For the 30,400 deg 2 T dwarf survey of Burgasser et 
al. (2003a) examined in Paper II, these simulations (using the stellar density normalization from Reid et al. 
1999) predict 22 and 45 disk T dwarfs with 700 < T e ff < 1300 for a = 0.5 and 1.5, respectively, independent 
of color constraints. These values straddle the current number count from this survey, roughly 36 T dwarfs 
(Burgasser, McElwain, & Kirkpatrick 2003; Burgasser et al. 2003a, 2004a, in prep; Tinney et al. in prep.). 

Surface densities for L and T dwarfs in the deep imaging case are substantially higher (^2 — 20 deg -2 
[100 K] _1 ) and approximately constant from 1000 < T e ff < 2500 K. This flattening is caused by the vertical 
extent of the Galactic disk, which truncates the surface density of warmer sources. For thinner disks (H z = 
200 pc), £ is smaller but increases toward cooler T e //. Hence, both the shape and magnitude of the surface 
density distribution can provide constraints on the vertical distribution of brown dwarfs for deep imaging 
surveys. Warmer M- and L-type halo stars and brown dwarfs, for which we assume a larger scale height (H z 
= 3 kpc), can rival their disk counterparts in surface density despite their lower space density. The coolest 
brown dwarfs (T e ff ^ 800 K) have a surface density distribution similar to the shallow survey case, as these 
objects are too dim to be detected beyond H z . For H z = 300 pc and a = 0.5, these simulations predict 2-3 
T dwarfs (500 < T eff < 1300 K) and 3-4 L dwarfs (1300 < T eff < 2300 K; Golimowski et al. 2004) in the 
160' 2 UDF z' field, sources that could be identified through followup deep imaging and/or spectroscopy. 



6. Summary 

Monte Carlo simulations of the substellar mass function have been presented, yielding LF and T e // 
distributions that can be directly compared to observations. A few salient points are worth reviewing: 

• Luminosity and T e // distributions for relatively simple realizations of the underlying mass function 
show a complex morphology in the brown dwarf regime, including: (1) a peak at the stellar/substcllar 
limit, (2) a paucity of sources in the L dwarf regime, (3) a rise in number densities for T-type and 
cooler brown dwarfs, and (4) a low-luminosity peak that depends on the minimum formation mass for 
brown dwarfs. 

• Variations in the stellar birthrate, minimum age, minimum formation mass, and choice of evolutionary 
model have minimal effect on the LF of T dwarfs, although L dwarf densities can be significantly 
skewed by any of these. As such, measuring the LF of the local T dwarf population provides the best 
means of constraining the substellar field MF. 

• Determining the minimum formation mass of brown dwarfs in the field will likely require the detection 
of significant numbers of objects with T e ff ^ 500, a temperature regime dominated by very low mass, 
non- fusing brown dwarfs. 

• Field L-type brown dwarfs, which evolve rapidly due to their higher luminosities, may be younger on 
average than field T dwarfs, a prediction that has some observational support (Gizis et al. 2000; Vrba 
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et al. 2004). For similar reasons, there may be a significant deficit of L-type dwarfs in the Galactic 
halo (r > 9 Gyr). 

• Unresolved multiplicity can enhance features in the observed LF, but these effects are generally small 
for binary fractions typical of brown dwarfs (10-20%). 

• Surface density estimates for the Hubble UDF suggest that a handful of L and T dwarfs will be present 
in that survey which can probe the disk scale height of brown dwarfs as well as detect a substantial 
number of halo low mass stars and brown dwarfs. 

The results presented here are qualitatively in agreement with those of Allen et al. (2004), who construct 
a two-dimensional grid of mass and age distributions to derive L and T e ff distributions for comparison 
via Bayesian analysis. In particular, many of the features in the LF identified in their simulations also 
appear here, despite differences in technique. Both studies therefore provide useful tools for constraining the 
substellar MF in the field. 

Paper II in this series will apply the simulations presented here to the local T dwarf LF derived from 
the 2MASS survey of Burgasser et al. (2003a), improving upon earlier estimates by Burgasser (2001) that 
were hindered by small number statistics. The simulations can also be used for a wide variety of imaging 
surveys, both as a predictive tool and as a means of probing the shape and scale of the substellar MF, the 
age distribution of cool halo dwarfs, the minimum "stellar" formation mass, and the vertical distribution of 
brown dwarfs in the Galaxy. 
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L. Moustakas during the preparation of the manuscript, and thanks the referee I. N. Reid for extensive 
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through Hubble Fellowship grant HST-HF-01 137.01 awarded by the Space Telescope Science Institute, which 
is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract 



A. An Analytic Correction for Space Density Measurements for Unresolved Multiple 



Unresolved multiple systems in a magnitude-limited survey can bias LF and space density measurements 
by hiding unseen members of the sample (underestimating source counts) and skewing photometry-based 
distance estimates for unresolved systems (overestimating space densities) . Since magnitude- limited surveys 
are commonly used to measure the LF, particularly in the field, an analytic expression for these effects is 
useful. 

Assume that the space density <p sys = J ^(Mbo^dMtoi has been measured for a large, shallow (i.e., 
ignoring disk scaleheight effects), magnitude-limited, and unresolved sample. Objects in this sample have 
an intrinsic binary fraction fbm and mass ratio distribution P(q), both of which are independent of mass, 
luminosity, age, etc. The space density can be represented as a sum over all (N) sources in the sample: 



NAS5-26555. 
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(Schmidt 1968), where 

Vmax — g a max — 3 V ' 

is the maximum volume sampled for a source with intrinsic brightness M to a limiting magnitude mu m 
over a surface area f2. Unresolved binaries in this sample with component brightnesses and F^ 2 \ and 
mass ratios q = /M^ , add additional sources to this sum, while also increasing the maximum distances 
(d max ) to which they can be detected. For each unresolved binary i, the correction to l/V max ,i oc d^ ax i is 

(l + ^Tyr 3/2 = (l + ?7)- 3/2 ^^ (A3) 
F i 

(see Eqn. 14), where it is assumed that L oc M 7 . The space density of individual sources, </w, can therefore 
be expressed in three terms: 

N sin N p „ N sec 

*»* = E^ + ^^ + ^^ (A4) 

*max.i ■ *max. i - v max,i 

% i i 

where N s i n = (1 — fbin)N is the number of single systems and N pr i — N sec = fbi n N are the number of 
binary primaries and secondaries, respectively. To the limit of large TV, variations in individual V max values 
can be averaged out and the summations replaced by Eqn. Al: 

4>md = (1 - fbin)4>sys + 2fbi n r(q)4>sys /»r\ 

= <f> sys (l-f bin (l-2T(q))) [ > 



where 

/•i 

P(q)(l + q~<)- 3/2 dq- (A6) 



r(«) = f 

Jo 



In the limiting case that all binaries have negligible mass secondaries (q — > 0), T(q) — > 1, and there is no 
correction to the distance estimates of the unresolved systems; all of the secondaries are added to the space 
density. For all binaries in equal-magnitude systems (P(q) = 5(q — 1)), T(q) = \ = 0.354, and 4> in d = 
(1 — O.293/bi n )0 S y S ; i.e., <f>i n d < 4> sys - This is due to the larger V ma x values for equal- magnitude systems, 
which overwhelms the increase in <\>ind from the addition of new secondaries. For the _P(^) distributions 
listed in Table 1 and diagrammed in Figure 13, and assuming 7 = 2.64 (Burrows et al. 2001), T(q) > 0.5, so 
that <j)i n d > 4> S ys- Table 5 lists values for <j>ind/4>sys for various combinations of fu n and P (<?)• These values 
are comparable to the simulated corrections diagrammed in Figure 14. 
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Tabic 1. Fundamental Distributions for Monte Carlo Simulations. 



Distribution Form 
(1) (2) 


Parameters 
(3) 


*(M) oc M~ a 


a = 0.0, 0.5, 1.0, 1.5, 2.0 


(logM-log M c ) 2 




oc e 2<r 2 


M c = O.1M ,<t = 0.627 a 


P(t) = b(T a -t) oc constant 




K e-( T <>-*)/r 8 


T = 10 Gyr, t 9 = 5 Gyr 


empirical 13 




((T -t)-t« ) 2 




<* Etl c 


N cl = 50, r cl = 10 Myr 


oc constant t < 1 Gyr 




P(Z) oc constant 


Z = Z Q 


P(q) oc constant 




oc e(9-!)/9c 


q c = 0.26 c 


from MF d 


a = 0.5 


a Parameters from Chabrier (2001). 

b Based on data from Rocha-Pinto et al. (2000). 

c Parameter fit from distribution of L and T dwarf binaries from Reid et 
al. (2001); Bouy et al. (2003); Burgasser et al. (2003b); and Gizis et al. 
(2003); see Figure 13. 

d Distribution based on Monte Carlo simulation of random pairings from 
an a = 0.5 MF; see also Kroupa & Burkert (2001). 
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Table 2. *(M) for Baseline Simulations' 1 . 



Mass (M ) 
(1) 


a = 




: 0.0 

2) 


a = 0.5 

(3) 


lognormal 
(4) 


a = 1.0 

(5) 


a = 1.5 

(6) 


a = 2.0 

(7) 


0.010-0.015 


0. 


.55 


1.5 


1.5 


4.3 


12 


33 


0.015-0.020 


0. 


.55 


1.3 


1.4 


3.0 


7.1 


17 


0.020-0.025 





.55 


1.1 


1.4 


2.4 


4.8 


9.9 


0.025-0.030 


0. 


.55 


1.0 


1.3 


1.9 


3.6 


6.6 


0.030-0.035 


0. 


.55 


0.94 


1.2 


1.6 


2.8 


4.7 


0.035-0.040 


0. 


55 


0.87 


1.1 


1.4 


2.2 


3.5 


0.040-0.045 


0. 


55 


0.83 


1.0 


1.2 


1.9 


2.8 


0.045-0.050 


0. 


.55 


0.78 


0.98 


1.1 


1.6 


2.2 


0.050-0.055 


0. 


55 


0.74 


0.91 


1.0 


1.3 


1.8 


0.055-0.060 





.55 


0.71 


0.83 


0.91 


1.2 


1.5 


0.060-0.065 


0. 


.55 


0.67 


0.81 


0.84 


1.0 


1.3 


0.065-0.070 


0. 


.55 


0.65 


0.74 


0.78 


0.92 


1.1 


0.070-0.075 


0. 


55 


0.62 


0.71 


0.73 


0.83 


0.92 


0.075-0.080 


0. 


.55 


0.61 


0.65 


0.68 


0.76 


0.82 


0.080-0.085 


0. 


.55 


0.59 


0.64 


0.64 


0.69 


0.73 


0.085-0.090 


0. 


.55 


0.58 


0.59 


0.60 


0.62 


0.65 


0.090-0.095 


0. 


.55 


0.56 


0.55 


0.57 


0.57 


0.59 


0.095-0.100 


0. 


55 


0.54 


0.53 


0.53 


0.53 


0.51 



a In units of pc 3 Mg 1 . Baseline simulations assume the Baraffe et al. (2003) models, 
0.01 < M < 0.1 Mq, 0.01 < t < 10 Gyr, and a constant birthrate. *(M) is normalized to 
0.55 pc — 3 Mg 1 averaged from the low-mass star MFs of Reid et al. (1999) and Chabrier 
(2001). 
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Tabic 3. $(M bo( ) for Baseline Simulations' 1 . 



1 \ 

Mbol (mag) 


oc — u.u 


a — u.o 


lognormal 


OL — l.U 


( V — ±.0 


OL — Z.U 




\ z ) 






(°) 


v°) 


(') 


9.0-9.5 


5.9o-7 


2.9c-7 


5.1e-7 


2.4e-7 








9.5-10.0 


1.5c-5 


1.7e-5 


1.6c-5 


1.8e-5 


1.9c-5 


1.4c-5 


10.0-10.5 


4.9e-5 


5.6e-5 


5.5e-5 


5.6e-5 


6.4e-5 


7.9c-5 


10.5-11.0 


0.00011 


0.00012 


0.00013 


0.00014 


0.00016 


0.00022 


11.0-11.5 


0.00021 


0.00023 


0.00027 


0.00030 


0.00034 


0.00044 


11.5-12.0 


0.00036 


0.00043 


0.00046 


0.00048 


0.00059 


0.00070 


12.0-12.5 


0.0018 


0.0019 


0.0019 


0.0020 


0.0021 


0.0024 


12.5-13.0 


0.011 


0.011 


0.011 


0.011 


0.012 


0.012 


13.0-13.5 


0.0083 


0.0088 


0.0094 


0.0096 


0.011 


0.012 


13.5-14.0 


0.0045 


0.0051 


0.0057 


0.0059 


0.0070 


0.0082 


14.0-14.5 


0.0042 


0.0048 


0.0054 


0.0057 


0.0070 


0.0088 


14.5-15.0 


0.0028 


0.0034 


0.0040 


0.0043 


0.0056 


0.0074 


15.0-15.5 


0.0029 


0.0036 


0.0042 


0.0047 


0.0062 


0.0089 


15.5-16.0 


0.0030 


0.0039 


0.0046 


0.0053 


0.0072 


0.011 


16.0-16.5 


0.0034 


0.0045 


0.0053 


0.0061 


0.0087 


0.013 


16.5-17.0 


0.0040 


0.0054 


0.0064 


0.0076 


0.011 


0.017 


17.0-17.5 


0.0053 


0.0071 


0.0086 


0.010 


0.015 


0.023 


17.5-18.0 


0.0061 


0.0086 


0.010 


0.013 


0.019 


0.030 


18.0-18.5 


0.0066 


0.0096 


0.012 


0.015 


0.023 


0.038 


18.5-19.0 


0.0062 


0.0098 


0.012 


0.016 


0.026 


0.045 


19.0-19.5 


0.0058 


0.0097 


0.012 


0.017 


0.030 


0.054 


19.5-20.0 


0.0052 


0.0094 


0.012 


0.017 


0.033 


0.065 


20.0-20.5 


0.0048 


0.0093 


0.011 


0.019 


0.039 


0.081 


20.5-21.0 


0.0041 


0.0087 


0.010 


0.019 


0.042 


0.094 


21.0-21.5 


0.0034 


0.0078 


0.0086 


0.018 


0.043 


0.10 


21.5-22.0 


0.0025 


0.0063 


0.0067 


0.016 


0.041 


0.11 


22.0-22.5 


0.0017 


0.0044 


0.0044 


0.012 


0.033 


0.091 


22.5-23.0 


0.00075 


0.0021 


0.0020 


0.0063 


0.018 


0.053 


23.0-23.5 


4.7o-5 


0.00014 


0.00013 


0.00044 


0.0013 


0.0040 



a In units of pc 3 mag 1 . Baseline simulations assume the Baraffe et al. (2003) models, 
0.01 < M < 0.1 M , 0.01 < t < 10 Gyr, and a constant birthrate. 



-24- 



Table 4. <f>(T e //) for Baseline Simulations' 1 . 



Teff (K) 


a = 0.0 


a = 0.5 


lognormal 


a = 1.0 


a = 1.5 


a = 2.0 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


onn ^nn 
zuu ouu 


n nnn9A 

U.UUUZ4: 


n nnn^Q 
u.uuuuy 


n nnnfiA 

U.UUUU4: 


n nn9i 

U.UUZl 


n nnfi9 
u.uuuz 


n m r 

U.UIO 


qnn Ann 


n nn9A 

U.UUZ4 


n nn<^9 
u.uuoz 


n nnfiA 

U.UUU4: 


n ni 7 

U.Ul I 


044 


12 


Ann— ^nn 

1UU ouu 


n nn^ 

u. uuoo 


0.0073 


n nns^ 

u. uuoo 


0.017 


0.039 


0.090 


^nn— finn 

ouu uuu 


0.0034 


0.0067 


0.0081 


0.013 


0.028 


0.057 


«nn 7nn 

UUU I uu 


n nn^ 
u.uuoo 


n nn^s 
u.uuoo 


n nn79 

U.UU I z 


n ni n 

U.U1U 


n ni q 
u.uj-y 


n n^7 
u.uo / 


7nn snn 

( uu ouu 


n nn^9 
u.uuoz 


n nn^9 
u.uuoz 


n nnfiA 

U.UUU4: 


n nn&7 

U.UUO I 


n ni ^ 

U.UFO 


n n97 

u.uz / 


eon nnn 
ouu yuu 


n nn^n 
u.uuou 


U.UU4:U 


n nn^fi 
u.uuou 


n nn7i 

U.UU I 1 


n ni 1 

U.U-L 1 


n mo 
u.uiy 


onn i nnn 
yuu— iuuu 


n nn9Q 
u.uuzy 


n nnAi 

U.UUOl 


n nnAQ 
u.uu^ty 


n nn^n 

U.UUOU 


n nnon 
u.uuyu 


n m a 


1 nnn 1 i nn 

IUUU - 1 1UU 


n nn9A 

U.UUZ^i 


n nn^^ 
u.uuoo 


n nnAn 
u.uu^tu 


n nnAS 

U.UU4o 


n nn^o 
u.uuuy 


n m i 

U.Ul 1 


1 1 nn 1 onn 

± 1UU — 1ZUU 


n nn9n 
u.uuzu 


n nn97 

U.UUZ ( 


n nn^^ 
u.uuoo 


n nn^s 

U.UUOO 


n nn^A 

U-UUD^ 


n nn&i 

U.UUOl 


i onn— 1 ^nn 

IZiUU -LOUU 


0.0015 


n nn9n 

U . UUZrU 


0.0024 


n nn9K 

U.UUZrO 


n nnAn 


n nn^Q 


1 ^nn 1 Ann 

lOUU I'lUU 


n nm 9 

U.UU-LZ 


n nm 7 

U.UUl I 


n nn9n 
u.uuzu 


n nn99 
u.uu^z 


n nn^i 

U.UUOl 


n nnAfi 

U.UU4U 


1 Ann 1 ^nn 

l^UU - lOUU 


n nn 1 1 

U.UUl 1 


n nni a 

U.UUl^t 


n nm & 
u.uuiu 


n nm q 
u.uuiy 


n nn9fi 

U.UUZO 


n nn^fi 

U.UUOU 


1500-1600 


0.00096 


0.0012 


0.0014 


0.0016 


0.0022 


0.0031 


1600-1700 


0.00088 


0.0011 


0.0013 


0.0015 


0.0019 


0.0026 


1700-1800 


0.00081 


0.00099 


0.0012 


0.0013 


0.0017 


0.0023 


1800-1900 


0.00077 


0.00095 


0.0011 


0.0012 


0.0015 


0.0021 


1900-2000 


0.00077 


0.00094 


0.0011 


0.0012 


0.0014 


0.0019 


2000-2100 


0.00097 


0.0011 


0.0013 


0.0014 


0.0017 


0.0021 


2100-2200 


0.0011 


0.0013 


0.0014 


0.0015 


0.0018 


0.0022 


2200-2300 


0.0011 


0.0012 


0.0014 


0.0014 


0.0017 


0.0020 


2300-2400 


0.0012 


0.0013 


0.0015 


0.0015 


0.0017 


0.0019 


2400-2500 


0.0025 


0.0027 


0.0028 


0.0029 


0.0032 


0.0035 


2500-2600 


0.0020 


0.0021 


0.0022 


0.0022 


0.0024 


0.0026 


2600-2700 


0.0031 


0.0032 


0.0032 


0.0033 


0.0034 


0.0035 


2700-2800 


0.0030 


0.0030 


0.0030 


0.0030 


0.0030 


0.0030 


2800-2900 


0.00022 


0.00024 


0.00025 


0.00025 


0.00028 


0.00030 


2900-3000 


0.00013 


0.00013 


0.00014 


0.00014 


0.00014 


0.00014 



a In units of pc 3 (100 K) . Baseline simulations assume the Baraffe et al. (2003) 
models, 0.01 < M < 0.1 M Q , 0.01 < t < 10 Gyr, and a constant birthrate. 



Table 5. Unresolved Multiplicity Correction Factor <f> md / '(j> svs for Various fbin and P(q) 



fbin 


P(q)<x5(q-1) 


1 


e (9-l)/9c 


MF Selection 


0.1 


0.97 


1.05 


1.01 


1.06 


0.2 


0.94 


1.10 


1.02 


1.13 


0.5 


0.85 


1.25 


1.05 


1.32 


1.0 


0.71 


1.51 


1.10 


1.64 
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(o) "f (b) [ (c) 




t (Gyr) t (Gyr) 



Fig. 1. — The five birth rates examined in this study (Table 1): (a) a constant birth rate; (b) an exponentially 
decreasing birth rate, with time constant r g — 5 Gyr; (c) a smoothed version of the empirical birth rate of 
Rocha-Pinto et al. (2000), with the burst events A, B, and C indicated (Majewski 1993); (d) a stochastic 
birth rate assuming star formation exclusively in N c i = 50 clusters randomly distributed over the age of the 
Galaxy, each described by a Gaussian birth rate distribution with half-width t c ; — 5 Myr; and (e) a halo 
birth rate in which all brown dwarfs are formed in the first 1 Gyr. The birth rates are related to the adopted 
age distributions by P{t) = b(T — t) = b(r), and are normalized such that 6(r)dr = T Q , where T a = 10 
Gyr. 
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Fig. 2. — Evolutionary tracks of T e ff versus time for the COND models of Baraffe et al. (2003, solid lines) 
and the clear atmosphere models of Burrows et al. (1997, dashed lines). Masses (from bottom to top, and 
labelled for the Baraffe tracks) of 0.001, 0.005, 0.01, 0.015, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.075, 0.08, 0.09, 
and 0.1 M are shown. Approximate locations for spectral types L0, L5, T5, ad T8 are indicated, based on 
empirical T e // determinations by Golimowski et al. (2004). 
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ognorma 
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Fig. 3. — Mass function distributions ('J'(M); number density per Solar mass) derived from the Monte 
Carlo simulations, for baseline parameters (P(t) = constant, 0.1 < t < 10 Gyr, and 0.01 < M < 0.1 M Q ). 
Distributions are sampled every 0.005 M Q . These results agree with the analytic forms to 5% or better, 
even for lower minimum mass limits. The normalization constant for the simulations, 0.55±0.18 pc -3 Mq 1 
at 0.095 M q , based on the low-mass star MFs of Rcid et al. (1999) and Chabrier (2001), is indicated by the 
solid circle. 
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o 0.010 
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Fig. 4. — Derived luminosity functions (&(Mboi); number density per magnitude) for the baseline MF 
simulations (Baraffe et al. (2003) models, 0.01 < M < 0.1 M Q , 0.01 < t < 10 Gyr, and constant birthrate). 
Distributions are sampled every 0.5 mag. The approximate location of spectral types LO, L5, T5, and T8 
are indicated, based on empirical M^oi determinations from Golimowski et al. (2004). 
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Fig. 5. — Derived T e // distributions (<J>(T e //); number density per 100 K) for the baseline MF simulations. 
Distributions are sampled every 100 K, and are slightly offset horizontally for clarity. The approximate 
location of spectral types LO, L5, T5, and T8 are indicated, based on empirical T e ff determinations from 
Golimowski et al. (2004). 



-30 - 



LO 



L5 



T5 



0.0100 



" 0.0010 



— All masses 

- - 0.075 < M < 0.100 Mo 
--- 0.012 < M < 0.075 M 
-- 0.001 < M < 0.012 M 



o.ooo- |-- 

3000 



2500 




Fig. 6. — T e // distribution for the a = 0.5 MF simulation with a lower-mass cutoff of 0.001 M Q , broken 
down into various mass bins: low mass stars (0.075 < M < 0.1 M Q ; dashed line), Deuterium-burning brown 
dwarfs (0.012 < M < 0.075 M Q ; dot-dashed line), and non-fusing brown dwarfs (0.001 < M < 0.012 M ; 
triple-dot-dashed line) . 
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Fig. 7. — Median mass versus Mm (left) and T e // (right) for baseline power-law MF simulations with a 
= 0.5 (black, open circles) and 1.5 (grey, filled circles). The uncertainty of the typical mass in each bin 
is indicated as the range of masses sampled by 63% of the population about the median, consistent with 
Gaussian ±lcr uncertainties. Data points are slightly offset horizontally for the a = 1.5 case for clarity. The 
Hydrogen and Deuterium burning mass limits (0.072 and 0.012 M , respectively) for the Baraffe models are 
delineated by dashed lines. 




Fig. 8. — Median age versus M^oi (left) and T e ff (right); symbols are those of Figure 8. Scatter in the ages 
for each bin are indicated and were calculated similarly to the mass uncertainties in Figure 8. Data points 
are slightly offset horizontally for the a = 1.5 case for clarity. 
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Fig. 9. — Comparison of &(Mb i) (left) and <j>(T e //) (right) for a = 0.5 (black lines) and 1.5 (grey lines) 
baseline simulations based on the Burrows et al. (1997, dashed lines) and Baraffe et al. (2003, solid lines) 
evolutionary models. Distributions are slightly offset horizontally for clarity. 




Fig. 10. — Comparison of $(M(, j) (left) and 3>(T e //) (right) for a = 0.5 baseline simulations for the 
five birth rates explored in this study. Distributions for each birthrate are slightly offset for clarity. The 
constant, empirical, and cluster birth rates show nearly identical distributions, while the exponential and 
halo distributions show significant variations in the Mboi ~ 15—17 {T e ff ~ 1200 — 2000 K; SpT T5-L3) 
trough and at faint luminosities/cold temperatures. Distributions are slightly offset horizontally for clarity. 
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Fig. 11. — Comparison of $(Mf, ;) (left) and $(T e //) (right) for a = 0.5 (black lines) and 1.5 (grey lines) 
baseline simulations for four minimum age limits, t — 0.001 (solid lines), 0.01 (dotted lines), 0.1 (dashed 
lines), and 1 (dot-dashed lines) Gyr. Distributions are slightly offset horizontally for clarity. 




Fig. 12. — Comparison of $(Mboi) (left) and $(T e //) (right) for a = 0.5 (black lines) and 1.5 (grey lines) 
baseline simulations for three minimum formation masses: M m i n = 0.001, 0.010, and 0.015 M©. Distributions 
are slightly offset horizontally for clarity. 
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q = M2/M1 

Fig. 13. — P(q) distributions employed to examine the effects of unresolved multiplicity: (dashed line) P(q) oc 
constant; (black line) P(q) oc e^ 9_1 ^ 9c ; (grey histogram line) P(q) from random pairing of primaries and 
secondaries drawn from an a = 0.5 MF with M min = 0.001 M Q . All distributions are arbitrarily normalized. 
The exponential distribution (with q c — 0.26) was chosen to fit an empirical mass ratio distribution (hatched 
histogram) constructed from 22 L and T dwarf binaries from Reid et al. (2001), Bouy et al. (2003), Burgasser 
et al. (2003b), and Gizis et al. (2003). 
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Fig. 14. — Multiplicity corrections ^md/^sys for a magnitude-limited survey that includes unresolved binary 
systems. Top panel: Variations as a function of binary fraction, fun — 0.1, 0.2, and 0.5, assuming a flat 
mass ratio distribution (P(q) = constant) and a = 0.5. Center panel: Variation between the three P(q) 
distributions employed (Table 1), assuming fun = 0.2 and a = 0.5. Bottom panel: Variations between two 
power-law MFs, a — 0.5 and 1.5, assuming P(q) = constant and fun — 0.5. 
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Fig. 15. — Surface densities (£; in units of deg~ 2 [100 K] _1 ) as a function of T e // for a shallow near-infrared 
survey limited to H — 16 Vega magnitudes (top) and a deep red optical survey limited to z' = 28 AB 
magnitudes (bottom). For the shallow survey, three cases are shown: a = 0.5 (thin black line) and a = 1.5 
(dashed line) disk populations, and an a = 1.5 halo population scaled by 0.3% (thick gray line). For the 
deep survey, disk scaleheight (H z ) values of 200 and 300 pc are employed for the disk populations and 3 kpc 
for the halo population. All models assume a lower mass limit of 0.001 M . 



